import numpy as np

def beam_element_stiffness_curved_3d(E, Ixx, Iyy, J, A, L_elem, R):
    # 计算局部刚度矩阵
    k_local = np.zeros((12, 12))

    EIxx = E * Ixx
    EIyy = E * Iyy
    GJ = E * J / (2 * (1 + 0.3))# 假设泊松比为0.3
    EA = E * A

    # 刚度矩阵的修正
    k_local[0, 0] = k_local[6, 6] = EA / L_elem
    k_local[0, 6] = k_local[6, 0] = -EA / L_elem

    k_local[1, 1] = k_local[7, 7] = 12 * EIyy / L_elem ** 3
    k_local[1, 5] = k_local[5, 1] = 6 * EIyy / L_elem ** 2
    k_local[1, 7] = k_local[7, 1] = -12 * EIyy / L_elem ** 3
    k_local[1, 11] = k_local[11, 1] = 6 * EIyy / L_elem ** 2

    k_local[2, 2] = k_local[8, 8] = 12 * EIxx / L_elem ** 3
    k_local[2, 4] = k_local[4, 2] = -6 * EIxx / L_elem ** 2
    k_local[2, 8] = k_local[8, 2] = -12 * EIxx / L_elem ** 3
    k_local[2, 10] = k_local[10, 2] = -6 * EIxx / L_elem ** 2

    k_local[3, 3] = k_local[9, 9] = GJ / L_elem
    k_local[3, 9] = k_local[9, 3] = -GJ / L_elem

    k_local[4, 4] = k_local[10, 10] = 4 * EIyy / L_elem
    k_local[4, 8] = k_local[8, 4] = 6 * EIxx / L_elem ** 2
    k_local[4, 10] = k_local[10, 4] = 2 * EIyy / L_elem

    k_local[5, 5] = k_local[11, 11] = 4 * EIxx / L_elem
    k_local[5, 11] = k_local[11, 5] = 2 * EIxx / L_elem

    k_local[5, 7] = k_local[7, 5] = -6 * EIxx / L_elem ** 2
    k_local[7, 11] = k_local[11, 7] = -6 * EIxx / L_elem ** 2
    k_local[8, 10] = k_local[10, 8] = 6 * EIxx / L_elem ** 2

    print("单刚", k_local)
    return k_local